/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package flowvis;

import flowsimulation.FlowSimulator;
import gui.IsoLinesRendererPanel;
import javax.media.opengl.GL;

/**
 *
 * @author waldo
 */
public class IsoLinesRenderer extends Renderer {
    
    // Original value taken from: Visualization.java
    public static float MAX_USER_DENSITY = 1f;
    public static float MIN_USER_DENSITy = 0f;
    private static int DIM;
    private static double[] rhos;
    private static int nRhos = 7;
    /**
     * 
     */
    private static boolean USE_SATURATION;
    /**
     * 
     */
    private static boolean USE_HUE;
    /**
     * Is the IsoLinesRenderer visible?
     */
    private static boolean VISIBLE = false;
    /**
     * Amount of levels to use in color banding
     */
    private static int NLEVELS = 127;
    /**
     * Scalar coloring.
     */
    static int scalar_col = 0;
    
    /**
     * GUI panel for this renderer.
     */
    private IsoLinesRendererPanel panel;
    
    /**
     * Constructor.
     * @param sim 
     */
    public IsoLinesRenderer( FlowSimulator sim ){
        super(sim);
        setVisible( false );
        useHue( false );
        useSaturation( false );
        panel = new IsoLinesRendererPanel( this );
        DIM = simulator.getGridSize();
        rhos = new double[nRhos];
        for( int n=0;n<nRhos;n++) rhos[n]=0;
    }

    /**
     * Assign zero value to the n+2 isolines.
     */
    public static void clearDensities(){
        for( int n=0;n<nRhos;n++) rhos[n]=0;
    }
    /**
     * Set the density value to the nth line
     * @param rho
     * @param n 
     */
    public static void setDensityN( double rho, int n ){
        if( n<0 || n>nRhos ) return;
        rhos[n] = rho;
    }
    /**
     * Get the density value from the nth line
     * @param rho
     * @param n 
     */
    public static double getDensityN( int n ){
        if( n<0 || n>nRhos ) return 0;
        else return rhos[n];
    }
    /**
     * @param new value for {@link IsoLinesRenderer#VISIBLE}.
     */
    public static void setVisible( boolean visible ){ VISIBLE = visible; }
    /**
     * @return value of {@link IsoLinesRenderer#VISIBLE}.
     */
    public static boolean getVisible(){ return VISIBLE; }
    /**
     * @return {@link IsoLinesRendererPanel} of this renderer.
     */
    public IsoLinesRendererPanel getPanel() { return panel; }
    /**
     * @return Current scalar coloring value.
     */
    public int getColoring() { return scalar_col; }
    /**
     * @param c , new <b>scalar_col</b> value.
     */
    public void setColoring(int c) { scalar_col = c; }
    /**
     * Sets the amount of bands (used to calculate color levels).
     * @param l , amount of levels.
     */
    public void setNLevels(int l){ NLEVELS = l; }
    /**
     * @return Amount of bands
     */
    public int getNLevels(){ return NLEVELS; }
    /**
     * @return 
     */
    public static boolean withHue(){ return USE_HUE; }
    /**
     * @param hue 
     */
    public static void useHue( boolean useHue ){ USE_HUE = useHue; }
    /**
     * @return 
     */
    public static boolean withSaturation(){ return USE_SATURATION; }
    /**
     * @param saturation 
     */
    public static void useSaturation( boolean useSaturation ){ USE_SATURATION = useSaturation; }
    /**
     * According to the scalar coloring value, the {@code rgb} values are calculated.<br>
     * 1. Switch to the operation which corresponds to the scalar value.<br>
     * 2.a {@link ColorMap#COLOR_BLACKWHITE}.<br>
     * 2.b {@link ColorMap#COLOR_BLACKWHITE_BANDS}.<br>
     * 2.c {@link ColorMap#COLOR_RAINBOW}.<br>
     * 2.d {@link ColorMap#COLOR_RAINBOW_BANDS}.<br>
     * 2.e {@link ColorMap#COLOR_BLACKBODY}.<br>
     * 2.f {@link ColorMap#COLOR_BLACKBODY_BANDS}.<br>
     * 2.g {@link ColorMap#COLOR_TEMPERATURE}.<br>
     * 2.h {@link ColorMap#COLOR_TEMPERATURE_BANDS}.<br>
     * 2.i {@link ColorMap#COLOR_GRADIENT_G2R}.<br>
     * 2.j {@link ColorMap#COLOR_GRADIENT_G2R_BANDS}.<br>
     * 3. Assign RGB color by calling {@code GL.glColor3f(R,G,B)}.
     * <br><br>
     * GL methods used:<br>
     * {@code glColor3f(red, green, blue)}: Specify new red, green, and blue
     * values for the current color.<br>
     * - GLfloat  red.<br>
     * - GLfloat  green.<br>
     * - GLfloat  blue.<br>
     * {@see http://www.opengl.org/sdk/docs/man/xhtml/glColor.xml}<br><br>
     * @param gl , related to OpenGL.
     * @param vy , float value to calculate the color to set.
     */
    private static void set_colormap(GL gl, float vy) {
        // Apply the color mapping equation
        if( vy>MAX_USER_DENSITY ) vy=MAX_USER_DENSITY;
        if( vy<MIN_USER_DENSITy ) vy=MIN_USER_DENSITy;
        vy=(vy-MIN_USER_DENSITy)/(MAX_USER_DENSITY-MIN_USER_DENSITy);
        float[] rgb = new float[3];
        switch( scalar_col ){
            case ColorMap.COLOR_BLACKWHITE:
                ColorMap.blackwhite( vy, rgb, USE_HUE, USE_SATURATION );
                break;
            case ColorMap.COLOR_BLACKWHITE_BANDS:
                vy *= NLEVELS;
                vy = (int) (vy);
                vy /= NLEVELS;
                ColorMap.blackwhite( vy, rgb, USE_HUE, USE_SATURATION );
                break;
            case ColorMap.COLOR_RAINBOW:
                ColorMap.rainbow( vy, rgb, USE_HUE, USE_SATURATION );
                break;
            case ColorMap.COLOR_RAINBOW_BANDS:
                vy *= NLEVELS;
                vy = (int) (vy);
                vy /= NLEVELS;
                ColorMap.rainbow( vy, rgb, USE_HUE, USE_SATURATION );
                break;
            case ColorMap.COLOR_BLACKBODY:
                ColorMap.blackbody( vy, rgb, USE_HUE, USE_SATURATION );
                break;
            case ColorMap.COLOR_BLACKBODY_BANDS:
                vy *= NLEVELS;
                vy = (int) (vy);
                vy /= NLEVELS;
                ColorMap.blackbody( vy, rgb, USE_HUE, USE_SATURATION );
                break;
            case ColorMap.COLOR_TEMPERATURE:
                ColorMap.temperature( vy, rgb, USE_HUE, USE_SATURATION );
                break;
            case ColorMap.COLOR_TEMPERATURE_BANDS:
                vy *= NLEVELS;
                vy = (int) (vy);
                vy /= NLEVELS;
                ColorMap.temperature( vy, rgb, USE_HUE, USE_SATURATION );
                break;
            case ColorMap.COLOR_GRADIENT_G2R:
                ColorMap.green2red( vy, rgb, USE_HUE, USE_SATURATION );
                break;
            case ColorMap.COLOR_GRADIENT_G2R_BANDS:
                vy *= NLEVELS;
                vy = (int) (vy);
                vy /= NLEVELS;
                ColorMap.green2red( vy, rgb, USE_HUE, USE_SATURATION );
                break;
            default: break;
        }
        gl.glColor3f(rgb[0], rgb[1], rgb[2]);
    }
    
    @Override public void visualize(GL gl) {
        if (!getVisible()) return;
        for( int n=0;n<nRhos;n++ )
          do_marching_squares( rhos[n], gl );
    }
    
    /**
     * for each cell ci of the dataset do
     *   int index=0
     *   for each vertex vj of ci do
     *     store inside/outside state of vj in bit j of index
     *   end for
     *   select intersection code from case table using index
     *   for all cell edges ej of selected case do
     *     intersect ej with isovalue v
     *   end for
     *   construct line segments from intersections
     * end for
     * @param isorho
     * @param gl 
     */
    private void do_marching_squares( double isorho, GL gl ){
        if( isorho<=0 ) return;
        //Total of cells, for gridsize=64: (64-1)^2 = 3969 
        boolean bv0, bv1, bv2, bv3;
        double v3, v2, v1, v0;
        double[] rho = simulator.getDensity();
        int mcase, i,j;
        double wn = getWinWidth() / (double) (DIM + 1);
        double hn = getWinHeight() / (double) (DIM + 1);
        double px0, px1, py0, py1, pxi, pyi, pxe, pye, t;
        for(j=0; j<DIM-1;j++) for( i=0; i<DIM-1; i++){
            mcase=-1;
            bv0=false; bv1=false; bv2=false; bv3=false;
            v3 = rho[i+j*DIM]; v2 = rho[i+j*DIM+1];
            v0 = rho[i+j*DIM+DIM]; v1 = rho[i+j*DIM+DIM+1];
            if( v3>isorho ) bv3=true; if( v2>isorho ) bv2=true;
            if( v0>isorho ) bv0=true; if( v1>isorho ) bv1=true;
            // Cases
            if( bv3==false&&bv2==false&&bv1==false&&bv0==false ) mcase=0;
            else if( bv3==false&&bv2==false&&bv1==false&&bv0==true ) mcase=1;
            else if( bv3==false&&bv2==false&&bv1==true&&bv0==false ) mcase=2;
            else if( bv3==false&&bv2==false&&bv1==true&&bv0==true ) mcase=3;
            else if( bv3==false&&bv2==true&&bv1==false&&bv0==false ) mcase=4;
            else if( bv3==false&&bv2==true&&bv1==false&&bv0==true ) mcase=5;
            else if( bv3==false&&bv2==true&&bv1==true&&bv0==false ) mcase=6;
            else if( bv3==false&&bv2==true&&bv1==true&&bv0==true ) mcase=7;
            else if( bv3==true&&bv2==false&&bv1==false&&bv0==false ) mcase=8;
            else if( bv3==true&&bv2==false&&bv1==false&&bv0==true ) mcase=9;
            else if( bv3==true&&bv2==false&&bv1==true&&bv0==false ) mcase=10;
            else if( bv3==true&&bv2==false&&bv1==true&&bv0==true ) mcase=11;
            else if( bv3==true&&bv2==true&&bv1==false&&bv0==false ) mcase=12;
            else if( bv3==true&&bv2==true&&bv1==false&&bv0==true ) mcase=13;
            else if( bv3==true&&bv2==true&&bv1==true&&bv0==false ) mcase=14;
            else if( bv3==true&&bv2==true&&bv1==true&&bv0==true ) mcase=15;
            
            switch( mcase ){
                case 0: // No intersection
                    break;
                case 1: // Check v0-v3, v0-v1
                    // Line v0-v3
                    px0 = wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0;
                    py1 = py0 - hn;
                    // Get intersection
                    t = (isorho-v0)/(v3-v0);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v0-v1
                    px0 = wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0 + wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v0)/(v1-v0);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 2: // Check v1-v0, v1-v2
                    // Line v1-v0
                    px0 = 2*wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0-wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v1)/(v0-v1);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v1-v2
                    px0 = 2*wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1=px0;
                    py1=py1-hn;
                    // Get intersection
                    t = (isorho-v1)/(v2-v1);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 3: // Check v0-v3, v1-v2
                    // Line v0-v3
                    px0 = wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0;
                    py1 = py0-hn;
                    // Get intersection
                    t = (isorho-v0)/(v3-v0);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v1-v2
                    px0 = 2*wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1=px0;
                    py1=py0-hn;
                    // Get intersection
                    t = (isorho-v1)/(v2-v1);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 4: // Check v2-v3, v2-v1
                    // Line v2-v3
                    px0 = 2*wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1 = px0-wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v2)/(v3-v2);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v2-v1
                    px0 = 2*wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1=px0;
                    py1=py1+hn;
                    // Get intersection
                    t = (isorho-v2)/(v1-v2);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 5: // Check v2-v3, v2-v1
                    // Line v2-v3
                    px0 = 2*wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1 = px0-wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v2)/(v3-v2);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v2-v1
                    px0 = 2*wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1=px0;
                    py1=py1+hn;
                    // Get intersection
                    t = (isorho-v2)/(v1-v2);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    // -------------Check v0-v3, v0-v1
                    // Line v0-v3
                    px0 = wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0;
                    py1 = py0 - hn;
                    // Get intersection
                    t = (isorho-v0)/(v3-v0);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v0-v1
                    px0 = wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0 + wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v0)/(v1-v0);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 6: // Check v2-v3, v1-v0
                    // Line v2-v3
                    px0 = 2*wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1 = px0-wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v2)/(v3-v2);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v1-v0
                    px0 = 2*wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0 - wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v1)/(v0-v1);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 7: // Check v3-v0, v3-v2
                    // Line v3-v0
                    px0 = wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1 = px0;
                    py1 = py0+hn;
                    // Get intersection
                    t = (isorho-v3)/(v0-v3);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v3-v2
                    px0 = wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1 = px0 + wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v3)/(v2-v3);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 8: // Check v3-v0, v3-v2
                    // Line v3-v0
                    px0 = wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1 = px0;
                    py1 = py0+hn;
                    // Get intersection
                    t = (isorho-v3)/(v0-v3);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v3-v2
                    px0 = wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1 = px0 + wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v3)/(v2-v3);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 9: // Check v2-v3, v1-v0
                    // Line v2-v3
                    px0 = 2*wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1 = px0-wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v2)/(v3-v2);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v1-v0
                    px0 = 2*wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0 - wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v1)/(v0-v1);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 10: // Check v3-v0, v3-v2
                    // Line v3-v0
                    px0 = wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1 = px0;
                    py1 = py0+hn;
                    // Get intersection
                    t = (isorho-v3)/(v0-v3);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v3-v2
                    px0 = wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1 = px0 + wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v3)/(v2-v3);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    //------------------- Check v1-v0, v1-v2
                    // Line v1-v0
                    px0 = 2*wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0-wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v1)/(v0-v1);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v1-v2
                    px0 = 2*wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1=px0;
                    py1=py1-hn;
                    // Get intersection
                    t = (isorho-v1)/(v2-v1);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 11: // Check v2-v3, v2-v1
                    // Line v2-v3
                    px0 = 2*wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1 = px0-wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v2)/(v3-v2);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v2-v1
                    px0 = 2*wn + (float) i * wn;
                    py0 = hn + (float) j * hn;
                    px1=px0;
                    py1=py1+hn;
                    // Get intersection
                    t = (isorho-v2)/(v1-v2);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 12: // Check v0-v3, v1-v2
                    // Line v0-v3
                    px0 = wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0;
                    py1 = py0-hn;
                    // Get intersection
                    t = (isorho-v0)/(v3-v0);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v1-v2
                    px0 = 2*wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1=px0;
                    py1=py0-hn;
                    // Get intersection
                    t = (isorho-v1)/(v2-v1);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 13: // Check v1-v0, v1-v2
                    // Line v1-v0
                    px0 = 2*wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0-wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v1)/(v0-v1);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v1-v2
                    px0 = 2*wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1=px0;
                    py1=py1-hn;
                    // Get intersection
                    t = (isorho-v1)/(v2-v1);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 14: // Check v0-v3, v0-v1
                    // Line v0-v3
                    px0 = wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0;
                    py1 = py0 - hn;
                    // Get intersection
                    t = (isorho-v0)/(v3-v0);
                    pxi = px0 + t*(px1-px0);
                    pyi = py0 + t*(py1-py0);
                    // Line v0-v1
                    px0 = wn + (float) i * wn;
                    py0 = 2*hn + (float) j * hn;
                    px1 = px0 + wn;
                    py1 = py0;
                    // Get intersection
                    t = (isorho-v0)/(v1-v0);
                    pxe = px0 + t*(px1-px0);
                    pye = py0 + t*(py1-py0);
                    // Draw the line
                    gl.glBegin(GL.GL_LINES);
                      set_colormap(gl, (float)isorho);
                      gl.glVertex2d(pxi,pyi);
                      gl.glVertex2d(pxe,pye);
                    gl.glEnd();
                    break;
                case 15: // All intersections
                    break;
            }
        }
    }
    
}
